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Escaping the accelerator; how, when and in what numbers do 
cosmic rays get out of supernova remnants? 
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ABSTRACT 

The escape of charged particles accelerated by diffusive shock acceleration from supernova 
remnants is shown to be a more complex process than normally appreciated. Using a box 
model it is shown that the high-energy end of the spectrum can exhibit spectral breaks even 
with no formal escape as a result of geometrical dilution and changing time-scales. It is 
pointed out that the bulk of the cosmic ray particles at lower energies must be produced 
and released in the late stages of the remnant's evolution whereas the high energy particles 
are produced early on; this may explain recent observations of slight compositional variations 
with energy. Escape resulting from ion-neutral friction in dense and partially ionized media 
is discussed briefly and some comments made on the use of so-called "free escape boundary 
conditions". Finally estimates are made of the total production spectrum integrated over the 
life of the remnant. 
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1 INTRODUCTION 

The idea that Supernova Remnants (SNRs) produce the bulk of 
the Galactic cosmic rays is now widely accepted. A key element 
in this general acceptance has been the development of a quite 
sophisticated theory of particle acceleration in shock waves, usu- 
ally called diffusive shock acceleration (DSA), which provides a 
compelling explanation of how a small number (of order 1CP 4 for 
conventional SNR parameters) of particles flowing into the shock 
can be accelerated to ultra-relativistic energies and become cosmic 
ray particles. Less attention has been paid to the question of how 
these particles then escape from the acceleration site and propagate 
into the interstellar medium although this is clearly an important 
question. However increasing attention is now being paid to this 
issue, in part because of theoretical developments relating to mag- 
netic field amplification which require it to be taken much more 
seriously, and in part because of interest in the idea that escaping 
particles illuminating nearby molecular clo uds might be sources 
of high energy gamma-rays and neutrinos (Casanova et al.|[2010l : 



Gabici & Aharonian 2007; Gabici et al. 2007 



2 THEORETICAL PRELIMINARIES 

In the standard simple theory of DSA particles are accelerated 
while within a diffusion length of the shock, there is no escape 
of particles upstream whatsoever, and the only escape downstream 
is by advection with the bulk flow jKrvmskiill 19771 ; lAxford et al.l 
ll977l : lBellll973 : lBlandford & Ostrikejl 19781) As long as the shock 
can be treated as an essentially planar structure propagating at fixed 
speed into an upstream medium capable of scattering particles this 



picture is certainly correct. Asymptotically a plane surface mov- 
ing at uniform speed, where the displacement is proportional to the 
time, will always overtake a randomly walking particle where the 
displacement grows only as the square root of the time. Thus at this 
level all particles are accumulated in the downstream region and 
there is no escape into the undisturbed far-upstream medium. 

The problem of course is that a lot of simplifying assumptions 
have been made in this picture. A real SNR shock is spherical and 
not planar; decelerating and not moving at constant speed; and the 
upstream scattering will certainly depend on the effective magnetic 
field strength as well as the wave spectrum and may even be totally 
suppressed in partially ionized regions of the ISM. The picture is 
further complicated when the nonlinear reaction of the accelerated 
particles on the shock structure, and the resulting secular evolution 
of the shock structure, is included. Of course, and relevant to this 
paper, if the scattering is totally suppressed particles can escape 
upstream even from a uniformly moving and infinitely extended 
planar shock, but it is important to note that this is a singular limit. 
Letting the shock radius go to infinity at finite scattering and let- 
ting the scattering go to zero at finite radius are two mathematical 
limiting processes that do not commute, so the zero scattering and 
infinite shock problem is formally ill-posed. For the problem of es- 
cape from supernova remnants what we actually want to study is 
the case of a large but finite shock radius and an upstream medium 
(the ISM in our Galaxy) which has a non-zero even if small level 
of scattering (as required by cosmic ray propagation models). 

Naively one may expect that, because 0.5 is larger than 0.4, 
a spherical blast wave expanding according to the Sedov scaling, 
R oc t 0,4 , between radius R and time t can be outrun by a randomly 
walking particle with displacement D growing as D oc t 0,5 so that 
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escape has to be considered a real possibility. In fact a more sophis- 
ticated argument is required because the particle in the course of its 
random walk may return to the shock even though on average it 
moves further away. It is in fact quite easy to show (see appendix) 
that a uniformly diffusing particle released at radius R\ from the 
origin will enter a sphere centred on the origin of radius Ro with 
probability Ro /R\ and escape to infinity without ever encounter- 
ing the sphere with probability 1 — Ro / R\ . Thus a particle located 
even as far upstream as ten shock radii will have a slightly larger 
than 10% chance of returning to the shock (slightly larger because 
the shock is expanding; it would be precisely 10% for a stationary 
spherical surface). This is of course a mathematical result which 
among other things depends on uniform and isotropic spatial diffu- 
sion in three spatial dimensions and as such must be treated with 
appropriate caution, but it does show that escape is not quite as 
straightforward as one might think. Further while escape is possi- 
ble in three dimensions it is formally impossible in one dimension; 
a particle randomly walking on a line visits every point infinitely 
often even though the recurrence intervals grow ever longer. Thus 
if the cosmic rays were effectively tied to a single magnetic field 
line they could never escape (at least formally). In reality of course 
there is some element of cross-field diffusion and the three dimen- 
sional result is probably more relevant. 

Even it there is a non-zero probability of return to the shock, 
for acceleration to continue this probability has to be very close to 
one (more precisely, it must be less than one by an amount of order 
U /c where U is the shock speed and c is the velocity of light). 
Thus there can be a large intermediate range where particles are 
no longer being accelerated, but have not really escaped from the 
shock either (if one defines this as having negligible probability of 
coming back to the shock front). This can be s hown more p r ecisel y 
using the so-called box model as was done in lDrurv et alj d2003h . 
For convenience we summarise the key points here while slightly 
generalising the treatment. 

2.1 An example of acceleration spectral breaks without 
escape in the "box" model 

The simplest approximate, but quite physically realisti c, treatment 
of sh ock acceleration is the so-called "box" model jDrurv et"aT1 
1999). In this the accelerated particles are assumed to be more or 
less uniformly distributed throughout a region extending one diffu- 
sion length each side of the shock, and to be accelerated upwards 
in momentum space at the shock itself with an acceleration flux 
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per unit surface area where U\ and L/j are the upstream and down- 
stream velocity and f(p) is the phase space density of the acceler- 
ated particles (assumed to have an almost isotropic distribution). 
If the diffusion length upstream is L\, and that downstream is L2, 
then 
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where m and K2 are the upstream and downstream diffusion coeffi- 
cients. To a first approximation we assume that both L\ and L2 are 
small relative to the radius of the shock and that we can neglect ef- 
fects of spherical geometry (in fact it is not too difficult to develop 
a spherical box model, but it unnecessarily complicates the argu- 
ment) so that the box volume is simply A(L\ + L2) where A is the 
surface area of the shock. The basic "box" model equation is then 
simply a conservation equation for the particles in the box; the rate 



at which the number in the box changes is given by the divergence 
of the acceleration flux in momentum space plus gains from injec- 
tion and advection and minus advective losses to the downstream 
region. 

^- t [A(L 1 + L 2 )4np 2 f(p)] + A^=AQ(p) 
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where Q(p) is a source function representing injection at the shock 
(only important at very low energies), Fx is a flux function rep- 
resenting advection of pre-existing particles into the system from 
upstream (normally neglected) and F2 is the flux of particles ad- 
verted out of the system and carried away downstream. The only 
complication we have to consider is that the box is time-dependent, 
with flow speeds, shock area and diffusion lengths all changing. 

The escaping flux is determined simply by the advection 
across the downstream edge of the box, that is 
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where we have to explicitly allow for the time varying size of the 
downstream region. Substituting this expression for F2 (p) and ne- 
glecting the advection of prexisting particles (the F\ (p) term) the 
box equation simplifies to: 
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3 dp Anp 2 

Partial differential equations of this form always reduce, by 
the method of characteristics, to the integration of two ordinary 
equations, one for the characteristic curve in the (p, t) plane 

dp _ Ui — U2 p 
di ~~ Li + L 2 3' 

and one for the variation of / along this curve 
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Q — everywhere except at the injection momentum and we can 
write the above equation as 

din/ _ 1_<M 1 <9Li Uj_ 

~ A dt 



(8) 



dt A dt L!+L 2 dt Lx + L 2 

But the shock area A is a function only of time so that 

dA _ A A 
dt ~ dt' 

and, although the upstream diffusion length does depend on both 
time and momentum, if we assume Bohm scaling for the two 
lengths so that 
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(where v is the particle velocity) we can write 
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< 1). Finally, noting that 
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and a fixed compression ratio, we can simplify equation l[8j to 

dln/ _ din A dhxtjhjh) SUx d lnp 
dt dt + dt Ux-U 2 dt '■ 

which integrates trivially to relate the value of / at the end of one 
of the characteristic curves, say at the point (pi, ti), to the value at 
the start, say at (to,Po)> as follows; 
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is the standard exponent of the steady-state power-law spectrum 
associated with shock acceleration. 

This rather beautiful result shows how the standard test parti- 
cle power-law is modified by a combination of effects as the box 
volume changes. As one would expect the amplitude varies in- 
versely as the shock area and also decreases if the upstream dif- 
fusion length (at fixed energy) increases, but with an exponent be- 
tween zero and one determined by the ratio of the upstream dif- 
fusion length to the total width of the diffusion region. It is very 
interesting that the result is not simply a variation inversely as the 
box volume, which one would naively expect from geometrical di- 
lution. This reflects the fundamental asymmetry between the up- 
stream and downstream regions, that upstream is empty outside the 
diffusion region whereas the entire downstream region is filled with 
accelerated particles. 

If we assume pure Bohm scaling the other differential equa- 
tion is also integrable so that the problem is reduced entirely to 
quadratures (of course only within the various approximations we 
are making; but still a remarkable result). Bohm scaling implies 
that the mean free path is of order and proportional to the particle 
gyroradius, so that if the particle charge is e 
pv 
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where a is a dimensionless parameter (of order ten). Substituting 
in the equation of the characteristic (equation (|6})we get 
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and noting the relativistic identity between kinetic energy T, mo 
mentum p and velocity v, 

- II 

dp ' 

we can integrate this as 
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(Ux - U 2 ) UxBx dt. 

For relativistic particles the kinetic energy and the momentum are 
essentially interchangeable with T — c^Jp 1 + m 2 c 2 —mc 2 « cp. 

Further progress can be made by assuming that the shock ve- 
locity follows the Sedov self-similar law for a strong spherical ex- 
plosion in a cold gas where the shock radius expands as R oc t 2//s 
and the shock velocity decreases as U oc t~ 3 ^ 5 . We further sup- 
pose that the effective magnetic field scales as a power of the shock 
velocity, 

B cB oc oc r 3M/5 . (21) 



The parameter /i characterises the degree of field amplification in 
the shock. No amplification corresponds to 11 = 0, energy eq uipar- 
tition would give fi = 1, and Bell's instability gives 
11 = 1.5. With these power-law scalings the integration is trivial 
and yields the acceleration paths, 
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These curves all rise extremely steeply, representing an initial 
phase of rapid acceleration, turn over, and then become asymptoti- 
cally flat. Physically it is clear that, as the shock slows and the field 
drops, the high energy particles cease to be significantly accelerated 
and simply diffuse further and further in front of the shock. 

A very important aspect of the curves is that they uniquely 
relate final energies (or equivalently momenta) to starting times. 
Asymptotically and at high energies the relation is a simple power- 
law; for Tx 2> To and to <C tx we have simply 
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Using this we can now translate the dilution factors from equation 
dl5b to additional power-law terms in the final momentum. Explic- 
itly, a given final momentum maps to a starting radius using a Sedov 
expansion-law: 

i?(to) oc ^ 5 oc p^™, (24) 
and thus the first term on the RHS of equation dl5t translates to a: 
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(note that we are thinking of the spectrum as a function of final mo- 
mentum pi at some fixed final time tx, so that A(tx) is constant). 
The final momentum can also be mapped to a starting velocity, us- 
ing again a Sedov expansion-law: 

[/(to) oct - 3/5 ocp? /(1+3fl) . (26) 
Thus the second term on the RHS of equation dl5| l scales as: 
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Ux(to)Bx(t ) J " \Ux(t ) 
Thus for the high-energy part of the spectrum we deduce that 
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Finally we need to determine the initial amplitude of / from 
an injection model. There are two main approaches to the injec- 
tion rate. The first, which is perhaps more consistent with the test 
particle approach, is to simply parametrise it by assuming that 
some fraction of the incoming thermal particles are "injected" as 
non-thermal particles at some suitably chosen "injection momen- 
tum" which separates the thermal particle population from the non- 
thermal. In other words one writes 



Q(p,t) = r](i)nxUx5(p - p lai (t)), 



(29) 



where nx is the upstream thermal particle number density, r\ <g 
1 is the injection fraction, pi n j is the injection momentum and as 
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usual 8 is Dirac's delta distribution. It should be clear that this is 
a parametrisation rather than a true injection model, however it, or 
equivalent parametrisations, have been very widely used, typically 
with rj taken to be a constant of order 10 -5 to 10 -4 for protons and 
Pinj ~ WrripUi where m p is the proton mass. However there is no 
real justification for this apart from the fact that it seems to yield 
reasonable results in many cases. 

With the above parametrisation the distribution function just 
above the injection energy can be simply determined by equating 
the acceleration flux to the injection flux, 



^Sl(U 1 -U2)f(p^)=vniU 1 , 
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On this injection model, noting that the injection momentum 
scales as the initial shock velocity, thus po oc U oc pj //< - 1+3 ' l ' ) , it 
is easy to see that the power-law exponent of the spectrum at high 
energies is modified to 
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dlnp Ui - U 2 1 + 3p, Ui - U 2 l + 3fi 

The first correction is essentially a hardening of the spectrum re- 
flecting increased injection at early times while the second is a soft- 
ening resulting from the combined dilution terms. 

The second approa ch adopts the idea, which can be traced 
back to the early work of Eichler ( 1979), that the injection process 
is inherently extremely efficient but that various feedback processes 
operate to reduce it to the point where the accelerated particles 
carry a significant part of the energy dissipated in the shock. For 
a standard spectrum close to p~ 4 the energy is almost uniformly 
distributed per logarithmic interval over the relativistic part of the 
spectrum. This suggests taking a reference momentum in the mildly 
relativistic region, po ~ mc, and determining / by a relation of the 
form 

^/(po)mc 2 a IpUi (f/i - U 2 ) (33) 



where A is a number which depends logarithmically on the upper 
cut-off and which for supernova remnants is probably somewhere 
between 10 and 100. 

On this injection model the reference momentum is fixed and 
the amplitude of the spectrum scales with the ram pressure of the 
shock, 
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It follows that the high-energy asymptotic slope is modified to 
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Again because of the high injection rate at early times the high- 
energy spectrum can even flatten, at least in theory. However it is 
important to note that this calculation is not strictly consistent in 
that we are assuming that the injection is sufficiently strong for non- 
linear effects to become important, but are ignoring the reaction 
effects on the shock structure in the acceleration calculation. The 
assumption that the diffusion length scales are small relative to the 
shock radius also breaks down and this will increase the importance 
of the dilution effects which are underestimated here. Thus both 
calculation should he viewed as toy models illustrating some of the 
effects that can occur rather than as definitive predictions. 



A further caveat is that both injection models are models 
for proton injection, the protons being the dynamically dominant 
species. Unfortunately very little is known about the factors con- 
trolling the injection of electrons and other minor species despite 
their importance for diagnostic tests. It has also been suggested 
that the injection may be nonuniform over the shock surface with 
a strong dependence on th e angle between the mean background 
field and the shock normal dVolk et alj2003l) . 

These results refer of course only to the asymptotic behaviour 
of the high energy part of the spectrum. As pi is decreased there 
comes a point where to is no longer small relative to ti. At this 
point all values of the final momentum map down to a small ap- 
proximately constant region and the spectrum becomes simply the 
standard spectrum with exponent —3Ui/(Ui — U 2 ). This break oc- 
curs at the point to which efficient acceleration is possible at that 
stage in the remnant evolution, and decreases as the the remnant 
ages. 

This calculation illustrates an important point. The break in 
the spectrum associated with the acceleration time scale becoming 
of order the dynamical time scale can be at a significantly lower 
energy than that at which particles are escaping from the shock and 
it is quite possible for a population of energetic particles to remain 
in the neighbourhood of the shock, and even return to it occasion- 
ally, without significant acceleration. This population will be geo- 
metrically diluted by the expansion of the shock and the increase 
in the upstream diffusion lengh (the region of space they fill con- 
tinues to expand faster than they can be supplied by the on-going 
acceleration) and this process can rather naturally produce spectral 
breaks at the point where the acceleration is no longer in equilib- 
rium with the expansion. The break in the spectrum reflects not just 
the geometrical dilution but also, and potentially of great interest, 
the time-history of injection at the shock. 

This effect also explains an apparent paradox in the simple 
theory. If the maximum particle momentum is calculated as 
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where i aC c is the acceleration time scale it is clearly a monotonic in- 
creasing function of time. However if one estimates the maximum 
momentum by arguing that the diffusion length has to be less than 
the shock radius, or equivalently that the local acceleration time has 
to be less that the dynamical time, this gives the condition 
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where R is the shock radius and R its expansion speed. Thus 
on this condition the maximum momentum is implicitly deter- 
mined by «:(p max ) « 0.1RR For a Sedov blast wave the product 
RR oc t~ ' 2 is a decreasing function of time (even if slowly) and 
as k is an increasing function of p one expects that p max should ac- 
tually decrease with time and not increase. The solution of course 
is that the spectrum has a break at the point given by equation d37b 
(where the acceleration is no longer able to be in equilibrium with 
the expansion), but continues to the higher energy given by the first 
condition of equation d36l l. 

To some extent it thus comes down to semantics. One can have 
a perfectly consistent picture of active SNRs being surrounded by a 
halo of particles which were accelerated at earlier times to high en- 
ergies, but which are no longer undergoing significant acceleration 
and are simply diffusing out into the Galaxy. But of course if this 
halo region were to become a large part of the Galaxy, and the prob- 
ability of further interaction with the shock vanishingly small, then 
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clearly it would make more sense to regard these simply as escaped 
particles. And of course once the shock becomes weak and the SNR 
effectively dies then all the particles, including those trapped in the 
interior, must eventually escape. Although very uncertain it is in- 
teresting to make some order of magnitude estimates. If the diffu- 
sion into the Galaxy occurs as indicated by studies of cosmic ray 
propagation, with a diffusion coefficient of order 10 28 cm 2 s _1 for 
particles around a GeV and an energy dependence something like 
E ' 5 , then PeV particles from a thousand year old SNR would dif- 
fuse a distance of order V3 x 10 41 cm « 180 pc. Thus particles 
produced in a hypothetical early pevatron phase of a historical SNR 
will not have escaped into the general Galaxy but will fill a roughly 
spherical region around the SNR of radius a few hundred parsecs. 
This is of course the motivation for looking for high-energy emis- 
sion from molecular cloud targets near SNRs. 



3 TRAPPING INSIDE THE SHOCK OF LOW-ENERGY 
PARTICLES 

In fact for the bulk of the cosmic rays, and in particular for those 
which dominate the energy density and for which the composition 
is well determined, trapping until the SNR dies is undoubtedly the 
better picture. At these relatively low energies the planar shock ap- 
proximation is good, the particles are accelerated at the shock and 
then swept downstream to accumulate in the interior of the rem- 
nant. Here of course they undergo adiabatic energy losses as the 
remnant expands, but the energy lost in this way goes to driving the 
shock and is thus recycled into the acceleration of new particles. 

The acceleration depends crucially on having strong scatter- 
ing in the neighbourhood (and in particular upstream) of the shock 
which produces a diffusion barrier at the shock preventing the es- 
cape of these low-energy particles. There are a wide variety of res- 
onant and non-resonant instabilities which can produce the neces- 
sary magnetic turbulence and it is generally believed that the lo- 
cal diffusion coefficient at the shock can be driven down to val- 
ues corresponding to Bohm scaling, that is a scattering mean free 
path of order the gyroradius. The sharp syn chrotron rims observed 
in young remnants, originall y in the radio i Achterberg et alj[l994l) 
and more recently in X-rays dYamazaki et alj|2004l) provide rather 
convincing direct observational evidence for such diffusion barriers 
associated with strong shocks. Were the diffusion typical of that in- 
ferred in the general ISM the rims would be much more extended 
and indeed as pointed out long ago shock acceleration would not 
be a viable mechanism for making the Galactic cosmic rays, e.g. 
iGinzburg & Ptuskinl fl98lh . Eventually however the shock will be- 
come too weak to sustain the required strong scattering, the diffu- 
sion barrier will collapse and at this point the particles in the in- 
terior will escape into the ambient medium as the shock (and by 
implication the SNR) dies. 



3.1 Energy-dependent composition changes 

This has the important consequence that at lower energies the cos- 
mic ray composition should be dominated by particles accelerated 
just before the remnant dies (thereby suppressing any freshly syn- 
thesised component c oming from the SN o r its progenitor wind). 
This was discussed in lDrurv&Keane]jl995h . The key point is that, 
as discussed above, the low energy particles trapped in the inte- 
rior are subject to adiabatic losses on the dynamical time scale of 
the remnant, and that this energy is then recycled into producing 



freshly accelerated low-energy particles. Thus the dominant low- 
energy particle population is that accelerated within the most re- 
cent dynamical time-scale of the remnant. This energy argument 
is the more fundamental one, but in fact even if one simply looks 
at the amount of matter swept up, far more matter is swept up at 
late times when the remnant is large than the amount (of order the 
ejecta mass) swept up prior to the onset of the Sedov phase. 

It is interesting to note that the picture is very different for 
the highest energy particles which have to be accelerated at the 
time when the shock is fastest (probably in the period leading up 
to the transition from free expansion to the Sedov phase) and it is 
tempting to speculate that this should produce slight compositional 
variations with energy. Indeed there appear to be such variations 
in recent experimental data, and we may well be seeing evidence 
of this effect. To put it succinctly, the composition at high energies 
should reflect acceleration early on prior to the mass sweep-up time 
when the ejecta interact with the same mass of swept-up material; 
the composition at low energies should reflect acceleration at the 
energy sweep-up time when the total thermal energy of the swept- 
up ambient material becomes comparable to the explosion energy 
(and the shock is weakening). In particular the apparently softer 
proton spectra compared to helium can be very easily explained by 
the older and weaker shocks propagating in a medium where more 
of the helium is neutral than in the strongly ionized environment of 
a young remnant. 



4 AN EMPIRICAL ESTIMATE OF MAGNETIC FIELD 
AMPLIFICATION 

As indicated in the introduction one reason that it is topical to dis- 
cuss these issues is because of the recent interest in magnetic field 
amplification in SNR shocks. Without field amplification the differ- 
ence between the two maximum momenta given by equations [36l 
and 137 1 above is quite small (at least for normal SNR parameters) 
and the discussion could be regarded as rather academic. The situ- 
ation is radically different if the effective magnetic field strength is 
a strongly increasing function of shock speed. In this case the local 
equilibrium energy becomes a strongly falling function of time and 
the gap between the two maximum energies can be many decades. 
Similar discussions can be found in ??. In fact if one wishes to ex- 
plain the particles between the "knee" in the cosmic ray spectrum 
at about 3 PeV and the "ankle" at 3 EeV in this way, then at least 
two decades have to be bridged in this way (chemistry can account 
for about one decade with the Iron spectrum extending to a higher 
total energy per particle than the proton spectrum). Thus observa- 
tionally we require that the product BRR be of order 10 17 V in the 
strong young shocks accelerating the highest energy particles and 
fall to a value more like 10 V just before the shock dies. 

This two decades decrease must be almost entirely in the mag- 
netic field strength, which must therefore be quite strongly depen- 
dent on the shock speed. Naively one may argue that the strongest 
shocks have speeds of order 10 4 kms -1 whereas the weakest are 
probably of order 10 2 km s 1 so that a roughly linear dependence 
fits well - but of course this is a very uncertain and crude estimate. 



5 SUPPRESSSION OF SCATTERING IN DENSE MEDIA 

A further complication is that the strong scattering required to 
maintain the diffusion barrier at the shock may fail in dense and 
partially ionized media because ion-neutral friction suppresses the 
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instabilities needed to disturb the magneti c field on the relevant 
scales ; a good recent discussion is that by Ptus kin & Zirakashvilil 
(120051). This was actually pointed out first in the seminal pap er by 
iBelll dl978h . Bell, and following him lDraine & McKed d 19931) as- 
sumed that at a finite distance in front of the shock the damping 
would totally suppress the wave scattering and that particles would 
then freely stream away at a fraction of the speed of light. Unfortu- 
nately both seriously overestimated the escaping flux by assuming 
that the intensity at the point where escape sets in would be the 
same as in the diffusion equation solution with scattering thr ough- 
out the upstream region. As pointed out in iDrurvet all (1996) one 
should actually solve the diffusion equation with a modified bound- 
ary condition of / = at the point where free streaming sets in 
which reduces the flux by a factor of order U/c. A similar transi- 
tion between a diffusion regime and a free streaming regime occurs 
in the theory of stell ar atmosphere s whe re this effect is well known. 

More recently iMalkov et al.l ( 120051) have advocated an inter- 
esting variant of this idea where instead of suppressing all the 
waves they point out that the damping will selectively suppress 
a range of wave-numbers which in turn produces an escape cone 
of finite opening angle in pitch-angle space. The authors consider 
the case of particles from a shock precursor beginning to enter a 
dense molecular cloud and argue that in the cloud there is a critical 
momentum p\ (their notation, and not to be confused with pi in 
the rest of this paper) such that particles with parallel momentum 
P|| | > pi are not scattered and escape at the speed of light. They 
assume that the effect of this is that only particles with Ipy | < p\ 
remain in the cloud to produce gamma-rays and that this automat- 
ically leads to a steepening of the spectrum by precisely unity be- 
cause the portion of phase space left is just pi/p. However if a 
significant part of the distribution function is removed, one can- 
not assume that the rest of the distribution function remains unaf- 
fected. In reality particles will rapidly pitch-angle diffuse into the 
loss-cone from parts of the distribution where the scattering is still 
strong and pull the distribution down towards zero. Malkov (per- 
sonal communication) has clarified that they regard the effect as a 
rapid and transient one, but it is far from clear that the time-scale 
for penetration of the precursor into the cloud is short relative to the 
pitch-angle relaxation time of the distribution. Nor is it clear what 
effect spatial variation of p\ and mirroring associated with spatial 
variation of the magnetic field strength have on the proposed mech- 
anism. 



6 A COMMENT ON FREE ESCAPE BOUNDARIES 

Although as argued above free-escape boundaries are expected to 
occur naturally in partially ionized media, they were originally in- 
troduced in discussions of shock acceleration for a totally different 
reason, namely as an artificial regularising procedure in the non- 
linear theory. The Monte-Carlo simulations of non-linear shocks 
pioneered by Ellison and his coworkers, at least in the versions de- 
veloped to date, require the modified shock structure to be planar 
and stat i onary . So also do the s emi-analytic theories developed by 
lEichled jl979t) , iMalkovl d 19971) and Blasi which depend crucially 
on the assumption of stationarity and planarity. As is well-known a 
steady modified planar shock produces such hard high-energy spec- 
tra that the accelerated particle pressure will diverge unless some 
cut-off is imposed (even a test-article f(p) oc p~ 4 spectrum will 
diverge logarithmically at high energies, and the non-linear effects 
make this much worse). In reality the location and shape of the cut- 
off is determined by factors such as the finite age and size of the 



shock, but in the calculations all of these are lumped into one artifi- 
cially imposed cut-off which is then adjusted to fit the system being 
studied. 

One way of imposing a cut-off is simply to assume that parti- 
cles vanish from the s ystem ("escape") when they reach some max- 
imum cut-off energy Ellison &Eichleill 19841) . This has the advan- 
tage of simplicity, but looks artificial and gives abrupt step-function 
cut-offs. An alternative which is widely used is to impose a so- 
called "free-escape" boundary at a finite distance upstream of the 
shock. The main advantage of this is that it gives a more natural- 
lookin g smooth cut-off and in addition, as argued bv lReville et al.l 
(2009), it is more physical in that the local diffusion barrier at the 
shock falls off at a finite distance upstream if there is strong lo- 
cal field amplification. This paper also gives a useful discussion 
and comparison of the two approaches. For a good recent discus- 
sion in the context of nonlinear shocks see Capr ioli et al.l d20ld) . 
A final point worth making is that steady models are by their na- 
ture incapable of including the geometrical dilution effects which 
we have seen to be important in determining the high-energy end 
of the spectrum. The dilution effects have to be approximated as 
escape, which leads to an over-estimate of the actual escape. 

While there is nothing wrong with the use of such artificial 
cut-offs, and indeed they are essential for the success of both the 
Monte-Carlo and semi-analytic models, their wide-spread use has 
perhaps contributed to impression that there is a strict dichotomy 
between particles that are being accelerated and particles that have 
escaped. The reality for high-energy particles, as argued above, is 
that there is no such sharp transition. Rather there is a gradual slow- 
ing of the acceleration and an increasing amount of time spent dif- 
fusing ever further upstream, at least as long as there is some finite 
level of scattering throughout the upstream region. It does of course 
depend on the level of upstream scattering and the scales that one 
is interested in. If the shock is small and the upstream scattering 
very weak, as is the case for particle acceleration at the Earth's 
bow shock for example, a description in terms of escape is entirely 
appropriate. If the scattering is almost entirely due to excitation of 
instabilities by the accelerated particles, and the pre-existing level 
of turbulence is very low, it is clear that a certain flux of escaping 
particles is needed to bootstrap the whole process (as observed in 
PIC simulations) and in this situation nonlinear effects could lead to 
a rather sharp transition in energy between trapped particles and an 
escaping population. Similar considerations apply to all numerical 
studies where the limited size of the computational box inevitably 
leads to the need to allow for escape. But as estimated above, for a 
SNR expanding into the ISM, the PeV particles even after a thou- 
sand years have probably only diffused a few hundred parsecs in 
the tangled ISM field whereas the SNR can have a radius of sev- 
eral parsecs to tens of parsecs; here it seems more appropriate to 
think of the accelerated particles as gradually falling out of equilib- 
rium with the shock acceleration and filling an extended halo than 
as abruptly escaping. 



7 THE TOTAL SOURCE SPECTRUM 

From the point of view of cosmic ray propagation theory what is 
actually required is an estimate of the source function for cosmic 
rays, in other words what is the total contribution per SNR to the 
galactic cosmic ray population? Leaving aside the complication of 
interacting remnants in superbubbles (although this must certainly 
occur, and indeed is indicated by the compositional data) we can try 
to answer this using the ideas developed above. A very similar dis- 
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cussion can be found in ICapriolietal] ( l2010h : see also lOhira et al.l 

boioh . 

Let us consider first the high-energy part of the spectrum 
(which we take to be the part between the 'knee' at about 3 x 
10 15 eV and the 'ankle' around 3 x 10 18 eV). The former, on the 
views presented here, is accelerated early on while the shock is 
fast and has a strongly amplified field. As the shock slows and the 
effective field decreases, the critical momentum, p„ at which the 
acceleration time-scale is comparable to the dynamical time-scale 
decreases and particles above this energy stop being accelerated 
and start to diffuse away from the remnant. For these particles the 
production rate is composed of two parts, the acceleration flux at 
the shock itself upwards through p* and the numbers of particles 
already accelerated which are left above p* as the cut-off momen- 
tum falls. These are respectively 



Qi 

Q-2 



f(p*)(Ui - U 2 )AtyR 2 



a 4nR J 
-47rp»/(p*)p,— — 



(38) 

(39) 
(40) 



where the second is an approximate estimate only (it assumes that 
the downstream particles uniformly fill the interior which is cer- 
tainly not the case, but it also neglects the contribution from the 
upstream population which should at least partially compensate for 
this). The ratio of the two terms is just 



'U 1 _U 2 
R 



(41) 



which, for the self-similar power-law evolution we are assuming 

oc t-d+Sfl/s 



with p, 

Qi = 2 
Q 2 l + 3n 



Ui 



, is just 

U 2 



Ui 



(42) 



Thus, as is clear on physical grounds, the first term dominates if 
there is no field amplification but the second term rapidly becomes 
more important as the field amplification increases. In general the 
two terms are of roughly equal importance. 

Particles at the critical momentum are produced for a period 
inversely proportional to the rate at which p* is decreasing, dt = 
—dpt/p*, and thus the spectrum released into the Galaxy is 



S(p*)(-dp,) = (Qi+Q 2 )dt- 



Qi + Q 2 



(43) 



Assuming power-law scalings so that the ratio of Qj to Q 2 is con- 
stant if follows that 



S(p*)oc-Q-ocplf(p*)R 3 . 
p* 

Further with these scalings 



R oc pi 



6/(1+3,*) 



(44) 



(45) 



and because acceleration from the injection energy po to p* is by 
definition fast so that there are no dilution effects, 



/(P.) = f(po) 



po 



-3U 1 /{U 1 -U 2 ) 



(46) 



If we now use the two injection models discussed above we deduce 
that, in the case of a constant injection number fraction, 



/( P0 )ocp - 3 oc P r^ 1+3 "> 



(47) 



dlnSjp* 
cHnp* 



= 2- 



3Ui 



+ 



Ui + 2U 2 \ 
1 + 3a* V Ui - U 2 J 



(48) 



which is excessively hard because of the very strong injection at 
early times. 

More interesting is the case of self-regulated injection where 
the scaling 



6/(l+3,i) 



f(Po) oc p 

exactly cancels the R 3 term and we deduce 



dhiSjp,) 
d In p, 



= 2 



3t/i 



Ui-Ui 



(49) 



(50) 



and thus the effective high-energy source spectrum has exponent 



In this case the escaping spectrum is the same as the production 
spectrum and energy i s equi-distributed ac ross both. A similar con- 
clusion is reached by Caprioli et alj J201(J) . 

In reality of course these toy models are far too crude and as 
noted above they should only be taken as an indication of some 
of the effects that can occur and as a guide to further numerical 
studies. They do however strongly suggest that the reason the spec- 
trum softens between the "knee" region and the "ankle" (if this is 
a source effect and not a propagation effect) is that the acceleration 
responsible for these particles occurs early in the remnant evolu- 
tion at a time when the shock is very fast, but not at full power. 
Prior to sweep-up the bulk of the explosion energy is in the form of 
kinetic energy and is not available for acceleration. Thus although 
it is possible to accelerate to very high energies, the total power is 
limited, and thus the source spectrum must turn down in this region. 
The smooth matching onto the standard spectrum at lower energies 
then follows naturally from the shock and remnant dynamics. 
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APPENDIX A: ESCAPE PROBABILITIES IN RANDOM 
WALKS 

Here we give the formal proof that a uniformly diffusing particle in 
three dimensions, if released at distance Ri from the origin, will es- 
cape to infinity without ever entering a sphere of radius Ro centred 
on the origin with probability 1 — Ro /Ri and will enter the sphere 
at least once with probability Ro/Ri- We use a standard technique 
in random walk theory and solve the steady state diffusion equa- 
tion for a source located at Ri and with an absorbing boundary 
condition on the sphere of radius Ro- The ratio of the flux being 
absorbed at Ro to that escaping to infinity is then the ratio of the 
respective probabilities as long as the random walk is well repre- 
sented by a diffusion process (this will be the case as long as the 
length scales are large compared to the scattering mean-free path 
A, that is J?i - Ro > A). 

The steady-state diffusion equation in spherical coordinates is 



dr 



( 2 df 



S(r-Ri) 



(Al) 



where as usual 8 denotes the Dirac delta distribution. We need to 
solve this with boundary conditions f(Ro) = and / — > as 
r — > oo. Away from the source the steady-state diffusion equation 
is trivially integrable to give 



f(r) 



(A2) 



where A and B are constants of integration. Considering the simple 
case where n is a constant this gives 



k r 



and thus the solution is, for Ro < r < Ri 
where C = A/k is another constant and 

f = ?± C (l M =1(^-A 

1 r \Ro Rj r \R 



(A3) 



(A4) 



(A5) 



for r > Ri. 

It follows immediately that the fluxes being absorbed at the in- 
ner boundary and escaping to infinity are in the ratio 1 to R\ / Ro — 1 
and thus that the probability of absorption at the inner boundary is 
i?o / Ri and of escape to infinity 1 — Ro /Ri . 

Note that the calculation can easily be generalised to the case 
where k is non-constant and space has n dimensions. In this case 
the solution inside the injection radius is just 



f(r) = A 



dr' 



l K(r') 



(A6) 



'Ro 

and that outside is 

f Hl dr' f°° dr' ( f 

f{r)=A j r'"-'«(r') X r'»-Mr') \J R 



dr' 



It follows that the probability of return to Ro if released at Ri 
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r dr' ( r 
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dr' 
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and of escape 
^ dr' 



'K(r') 



f Rl dr' ( r dr' \ 
J Ro r'"-*K(r>) \J Ro r"-'K(r')J 
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As is intuitively clear escape is easier the more spatial dimen- 
sions are available and the more k is an increasing function of ra- 
dius, but it is interesting that this can be quantified in such simple 
closed formulae. 



